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Abstract 
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We determine the elastic moduli of two-dimensional assemblies of 
^ , disks by computer simulations. The disks interact through elastic con- 

Q>^ I tact forces, that oppose the relative displacement at the contact points 

OO ■ by means of a normal and a tangential stiffness, both taken constant. 

Our simulations confirm that the uniform strain assumption results in 
inaccurate predictions of the elastic moduli, since large fluctuations in 
f^^ ■ particle displacements and rotations occur. We phrase their contribu- 

^^ . tion in terms of the relative displacement they induce at the contact 

points. We show that the fluctuations that determine the equivalent 

continuum behavior depend on the average geometry of the assembly. 

We further separate the contributions from the center displacement 

5^ I and the particle rotation. The fluctuations result in a relaxation of the 

system, but along the tangential direction the relaxation is generally 
entirely due to rotations. We consider two theoretical formulations 
for predicting the elastic moduli that include the fluctuations, namely 
the "pair-fluctuation" and the "particle-fluctuation" method. They 
are both based on the equilibrium of a small subassembly, which is 
considered representative of the average structure. We investigate the 
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corresponding predictions of the elastic moduli over a range of coordi- 
nation numbers and of ratios between tangential and normal stiffness. 
We find a significant improvement with respect to the uniform strain 
theory. Furthermore, the dependence of the fluctuations on coordina- 
tion number and ratio of tangential to normal stiffness is qualitatively 
captured. 

Keywords: granular media, equivalent continuum, elastic moduli, fluctua- 
tions, DEM simulations 

1 Introduction 

This study investigates the elastic mechanical properties of dense random 
isotropic two-dimensional assemblies of disks. Our work is framed in the con- 
text of micromechanics, which focuses on the relation between macroscopic 
behavior and microscopic interactions. At the macroscopic scale, stress and 
strain are measured, whose relation is determined by the elastic moduli of 
the equivalent continuum. At the microscale, forces arise between contacting 
particles that for quasi-static deformations must satisfy the balance of force 
and moment for all grains. 

In the context of elasticity, contact forces oppose the relative displace- 
ment between contacting grains by means of a contact stiffness. The work of 
Poritzky [Ij about contacting thin disks reveals for its normal component a 
dependence on the overlap that rapidly weakens with the confining pressure. 
We assume for our assemblies the limit case of constant normal stiffness. 
Since we focus on small displacements, we also consider the tangential stiff- 
ness constant [21 E]- By ignoring all inhomogeneity in the contact stiffness, 
we specifically focus on the effects of geometric disorder. 

Given the geometry of the contact network, the constitutive law for the 
contact forces determines the macroscopic mechanical properties once the rel- 
ative displacements at the contact points are known. Many micromechanical 
studies in El m El [HI El [ini El [12] use the so-called 'average strain hypoth- 
esis', where the relative displacements reduce to the contribution from the 
average strain. The resulting prediction of the elastic moduli is inaccurate 
[T3l [lH [15] , particularly for the shear modulus, as in general grains undergo 
additional displacements in order to attain equilibrium. More sophisticated 
predictions have recently been developed that incorporate such fluctuations, 
in [TB] for three-dimensional systems and in [T7] and [IB] for two-dimensional 



ones. They are based on the idea that even though the deformation is not 
uniform, its fluctuations are strongly correlated. The study by Koenders [2U] 
suggests that they occur with correlation lengths in the order of a few di- 
ameters, which implies that subassemblies of such a size already contain the 
essential features of the global structure. 

The validity of the theoretical predictions in [T7] and [TS] is investigated 
here by comparison with the elastic moduli computed by means of Discrete 
Element Method (DEM) simulations. The analysis is performed for various 
ratios of tangential to normal stiffness and coordination numbers. We con- 
sider assemblies with both larger and smaller coordination number than the 
onset of iso-staticity for disordered frictionless systems. This onset equals 
4 in two dimensions and is recurrent in numerical simulations, as dense as- 
semblies are usually obtained by means of an initial frictionless compression. 
In experiments |T9] on hard disks, coordination numbers smaller and larger 
than 4 are due, in turn, to the presence of friction and of ordered structures. 
In DEM simulations, disordered assemblies with larger coordination numbers 
than the onset of isostacity can be obtained, for example by neglecting in 
the constitutive law for the contact force the increase in the normal stiffness 
with the interpenetration or by applying a large pressure. Such systems are 
mainly of academic interest, particularly in the development of statistical 
approaches aimed at predicting the evolution of disordered systems. Also in 
three dimensions the scientific literature concentrates on samples with larger 
coordination numbers than the frictionless onset of iso-staticity, which equals 
6 in that case. However, studies concerned with the issue of numerically re- 
producing experimental results |22] , [23] emphasize that lower ones might be 
relevant to practical purposes. 

The outline of the study is as follows. Firstly, the basic micromechanical 
quantities of interest are defined in Section [21 Then, a concise description is 
given in Section [3] of theoretical approaches for predicting the elastic moduli 
based, in turn, on the average strain assumption and on the inclusion of 
displacement fluctuations. The performed DEM simulations are described in 
Section m The corresponding results are analyzed in Section O and compared 
to the results the theoretical predictions. The final section is dedicated to 
discussion of the results. 



2 Micromechanics 

We consider the contact between disks p and q of radius BF and -R^, respec- 
tively. The contact is identified by the unit vector nf' that points outwards 
from p along the line that joins the centers. The unit vector tf' is tangent to 
the contact (see Figure [1]). In components, 

nP« = (cos^P^sin^P^), 
tP" = (-sin^P^cos^P"), 

where 9^"^ is the contact orientation, counted counterclockwise from the hor- 
izontal axis. For future reference, we also define the branch vector /f , that 
joins the centre of particle p to that of particle q pointing outwards from p, 
i.e. 

/f = {RP + Ri) nf . 

As in monodisperse two-dimensional assemblies crystallization occurs, a log- 
normal distribution for the particle radii is adopted. Contacting particles 




Figure 1: Contact geometry: contact orientation and normal and tangential 
vectors to the contact. 

interact by means of contact forces. We denote by /f the i-th component 
of the force exerted on particle p by particle q. It has normal and tangential 
component to the contact /^'^ and /f, i.e. 



pq — fPQ^PI 
n ~ Ji "-i ' 

PQ _ fPQ+PI 



rpq _ rpq^pi 



In the hypothesis of quasi-static deformations, contact forces satisfy the bal- 
ance of force and moment on each grain. For instance, on particle p, 
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0, 

where the sum is over particles q that are in contact with it and Cjk is the 
two-dimensional permutation tensor. Contact forces oppose the relative dis- 
placement between contacting particles by means of a contact stiffness, whose 
normal and tangential component we denote, in turn, by kn and kt, both 
considered constant. If we denote by Af^ the relative displacement between 
particles p and q that are in contact, and by A^^ and Aj"' its normal and 
tangential component, we have 

fpq — I. /\pq 

Jn ~ '^n'-^n i 

/r = kAi\ 

/r = /fnf + /rc- (1) 

The theory of contact elasticity predicts the decay with the distance from 
the contact zone of the effects of contact interactions. For small enough 
deformations and stiff enough particles, as we assume to be the case here, 
contact interactions can be assumed to be confined to a contact point, and 
the grains kinematics to be approximated by that of rigid bodies. Therefore, 



Af = [/f - f/f + dj (i?^fi« + RPQF) n; 



pq 



where Uf and U^ are the displacement of the centre of particle p and g, VP 
and VL'^ their rotation. 

At the macroscopic, continuum level, the relevant quantities are the stress 
tensor aij and the strain tensor e^j, whose components are taken positive in 
compression. Contact forces and the geometry of the particle arrangement 
determine the expression of the former, as [2] 



^"• = ^E E m^ (2) 



09 ceciOg) 



that is the average over the area of interest 5* of Cauchy's stress [27]. The 
orientation 9g varies between and vr. Finally, the superscript pq referring 



to contacting particles has been replaced by for corresponding contacts c. 
Expression ([2]) emphasizes the dependence of the macroscopic behavior on 
the average force over equally oriented contacts. Experimental observations 
and numerical simulations [281 ISH ED] suggest for it the same dependence 
on the contact orientation as for the effects of the average strain. That is, 
compressive forces are the larger the closer the contact orientation is to the 
direction of major compression, and tangential forces have their maximum 
at contacts oriented at 45 degrees from it. 

Primary geometrical characteristics of granular assemblies are coordina- 
tion number F, i.e. the average number of contacts per particle, contact 
density ^5, i.e. the average number of particles per unit surface, and the 
contact distribution function E{9) [HI], defined such that E{6)d6 gives the 
probability of finding a contact with orientation 9 in the interval (6', 9 + d9), 
9 G (0, tt). In the case of isotropic assemblies, as considered here, the contact 
distribution function becomes E{9) = l/vr. With the use of such quanti- 
ties and by denoting averages over equally oriented contacts by overbars, 
expression (^ transforms into the integral 
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r Ei9)fM0)d0. (3) 

^0 



The elastic stiffness tensor relates stress and strain. In isotropic systems, it 
is fully described by the effective bulk modulus K and shear modulus G, as 

CTu + 0-22 = 2K(eii + £22) 

CTll - Cr22 = 2G(eii - 622). (4) 

We determine them by means of DEM computer simulations, by applying, 
in order, an isotropic compressive deformation ef, and a shear deformation 






.K_l I 0\ n MO 



e^i - eo ( Q M , e^j - eo ( Q _M , (5) 

and measuring the corresponding stress response. By eo we denote a magni- 
tude of the imposed strain. 



3 Theoretical modeling 

If the deformation e^j is prescribed, the theoretical prediction of the elastic 
moduli requires that of the corresponding stress tensor aij. Due to ([T]), this 



in turn requires a kinematic localization assumption, that is, the expression 
of the relative displacements between contacting particles as a function of 
eij. This procedure is depicted in Figure [2l Various kinematic localization 
assumptions are considered in this section, such as the uniform strain as- 
sumption and more sophisticated approaches that account for fluctuations. 



Macroscopic 
level 



Microscopic 
level 



Stress 



Strain 



Averaging 



Kinematic localisation 
assumption 



Force 



Contact constitutive 
relation 



Relative 
displacement 



Figure 2: Kinematic localisation assumption. 



3.1 Uniform strain 

In the hypothesis of uniform strain the relative displacement between con- 
tacting particles becomes 

^pg ^ ipg 

Its average over equally oriented contacts in the case of isotropic compression 
is isotropic. Only its normal component A^'*'^ differs from zero, namely 



eo^ 



(6) 



where I is the average length of the branch vector. In the case of shear, the 
average normal and tangential component A*^ '^"'^ ^"^ ^^ ^^^ 
oriented contacts are, in turn. 



and Aj of A^^ over equally 






e()lcos29, 
—eolsm29. 



(7) 



For isotropic assemblies of disks, the average strain assumption results in the 
bulk and shear moduli K"^ and C^ |0: 



K' _ nsT 

G' _ nsT 
kn 16 
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P, (8) 



that are upper bounds. The symbol P denotes the average over all contacts of 
the squared length of the branch vector. As discussed in [32], polydispersity 
makes the average length of the branch vector differ in general from the 
average diameter. For our size distribution I is about 3% larger than the 
average diameter, and P is about 15% larger than the average diameter 
squared. 

3.2 Approaches that incorporate displacement fluctu- 
ations 

Two approaches are considered that incorporate the fluctuations into the 
prediction of the elastic moduli, namely the particle-fluctuation (IPF) and 
the pair- fluctuation method (PF). They are discussed in detail in [T7] and in 
|18j . respectively. In [17], the IPF approach is applied numerically and proven 
to give upper bounds to the effective moduli, while in [18j the analytical 
solution that corresponds to the PF method is presented. Both approaches 
deal with the issue of equilibrium of a small assembly, made of a chosen 
particle A and of a contacting pair AB, respectively, surrounded by their 
first neighbors. Therefore, in these models fluctuations are assumed to occur 
with correlation lengths in the order of 3 or 4 diameters. 

In the more general case, the relative displacement between contacting 
particles can be decomposed into the contribution from averages and fluctu- 
ations, that is, from the imposed strain and the average particle rotation, on 
one side, and from fluctuations in particle displacements and rotations, on 
the other. In the absence of macroscopic rotation, as it is the case here, in 
an isotropic system the average particle rotation is zero [33]. If we denote 
by uf and u'^ the fluctuation of contacting particles p and q in the center 
displacement, and by uj^ and u'^ the fluctuations in rotations, we can write: 




Figure 3: Small assembly considered in the "particle- fluctuation" method, 
centered on particle A. The deformation of the particles shown in gray is 
according to uniform strain, while A is allowed to fluctuate. Its fluctuations 
stem from the balance of force and moment on it. 

In the IPF and the PF approach, only the chosen particle or pair is 
allowed to fluctuate, while their neighborhood is compelled to move according 
to the average strain. This is depicted in Figure [3] for the IPF case. As a 
result, in the IPF method the relative displacement between particle A and 
its r-th neighbor reduces to 

/\Ar _ lAr ~A tdA-A.At 

SO that the three equilibrium equations for particle A can be solved for its 
three unknown fluctuations. In the pair-fluctuation method, 

\AB lAB I /~B ~A\ (ryA-A , ryB ~ B\ ,AB 

but 

Af = e,,lf^ - uf - R^uHf , 
Af = 6,,Zf - nf - i?^i;^^^ 

where r denotes a neighbor of particle A different from i?, and s a neighbor 
of particle B different from A. The six equilibrium equations of particles A 



and B can be solved for the corresponding six fluctuations. In both cases, 
the small size of the solving system permits an analytical formulation. 

In this work, the equilibrium equations are solved numerically, particle 
by particle and pair by pair for the IPF and the PF approach, respectively, 
for the corresponding fluctuations. Using the procedure sketched in Figure 
121 the corresponding contact forces are computed and the stress and elastic 
moduli estimated. The resulting fluctuations are also analyzed in terms of 
the associated deformation mechanisms in the sections that follow. 

4 DEM simulations 

DEM simulations have been performed of large isotropic assemblies of 50,000 
disks, with radii from a lognormal distribution. Periodic boundaries have 
been employed to reduce boundary effects. 

In the DEM method, the deformation of the assemblies is computed by 
numerically integrating in time the equations of motion of all particles, and 
the results can be employed to analyze the actual deformation mechanisms. 

Initial equilibrium states with coordination numbers F = 3.5, F = 4 and 
F = 5 have been first prepared. To obtain the first one, a succession of 
frictionless and frictional compressions has been employed. The second one 
results from the isotropic compression of a frictionless gas, and the third one 
results from an additional isotropic compression. 

The assemblies obtained this way have been subjected to the two strain 
paths specified by eqn.([5]). Then, the corresponding stress response has been 
computed in order to determine the effective bulk and shear moduli K and G. 
Such simulations have been performed with bonded contacts: that is, neither 
contact creation nor disruption has been considered, which is appropriate for 
studying elastic behavior at small strains. At each coordination number, 
different ratios h/kn have been used, in the range between 0.05 and 1.0. 

5 Micromechanical analysis 

In this section, the results from the DEM simulations and the theoretical ap- 
proaches introduced in Section [3] are analyzed. The comparison between the 
elastic moduli from the DEM simulations and the uniform strain assumption 
emphasizes the relevance of displacement fluctuations. The quality of the 
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estimates of the elastic moduli given by the IPF and the PF methods is also 
evaluated. Then, the deformation mechanisms induced by the fluctuations 
at the macroscopic scale are analyzed, and the performance of the IFF and 
PF approaches interpreted in terms of their capability of capturing them. 

5.1 Moduli 
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Figure 4: Bulk modulus K: from DEM simulations (solid line), uniform 
strain assumption (dotted line), IFF (dashed) and PF (dashed-dotted) ap- 
proaches. Results for three coordination numbers F: F = 3.5 (D), F = 4 (A), 
F = 5 (O) aiid various stiffness ratios kt/kn- Bottom right: estimate perfor- 
mance for uniform strain assumption (dotted) and IFF approach (dashed), 
for F = 3.5, 4.0, 5.0. 

The bulk and shear moduli that result from the DEM simulations and 
the theoretical approaches discussed in Section [3] are compared in Figures 
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HlandO They are presented in the dimensionless form K/k^ and G/kn, as 
function of the ratio kf/kn and for the three chosen coordination numbers. 
In the fourth frame of the two figures, the quahty of the estimates is plotted 
in terms of their ratio to the effective moduh. The average strain prediction 
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Figure 5: Shear modulus G. Same symbols as in Figure HI 

performs poorly, especially at small coordination numbers. The inclusion 
of the fluctuations with as few degrees of freedom as in the IFF and PF 
methods does not guarantee an improvement with respect to the average 
strain assumption, as the case of the bulk modulus at F = 5 proves. 

The DEM simulations point out the dependence of the bulk modulus 
on the tangential stiffness, ignored by the average strain assumption, but 
captured once fluctuations are accounted for. The reason tofor such a de- 
pendence lies in the fact that even though the tangential forces are zero on 
average in hydrostatic compression, they are required at the particle level 
in order for the forces and moments to balance. By opposing the relative 
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displacements along the tangential direction, the tangential forces stiffen the 
assembly, the more so the larger the tangential stiffness. The general belief 
that the average strain assumption gives a good approximation of the bulk 
modulus seems hardly acceptable, at least below the onset of iso-staticity for 
frictionless systems and at small ratios kt/kn- Regarding the shear modulus, 
the inclusion of the fluctuations into the prediction results in a significant im- 
provement with respect to the average strain assumption, although important 
deviations from the DEM simulations remain, especially for low coordination 
number and low stiffness ratios kt/kn- 

5.2 Relative displacements 

At the macroscopic level, only the average over equally oriented contacts 
of the relative displacement between contacting particles is important, as 
follows from ([1]) and ([2]). We label by Af' the relative displacement between 
contacting grains p and q due to their fluctuations: 



and by A^ and A^*^ its normal and tangential component. 



pq 






"^m 



When an isotropic compression is applied, only the former contributes to 
the stress. Its distribution is uniform and has average A„, such that 

A„ = f3nK. (9) 

where A^ is given by expression (El). As the fluctuations relax the system, 
I3n is negative. It follows after some algebra from ([2]), the first of ([8]) and ([9]) 
that 

K = {l + l3n)K\ (10) 

Corresponding values of K'^ / K are shown in the fourth frame of Figure (jl]). 
Typical group averages for the case of shear loading are shown in Figure 
ini where along the tangential direction the contributions from the center 
displacement and the rotation have been separated. The average normal and 
tangential relative displacements induced by the fluctuations over equally 
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Figure 6: Relative displacements between contacting particles averaged over 
equally oriented contacts: normal component (D), tangential component due 
to center displacements (A) and rotations (x). Results for shear deformation 
at coordination number F = 3.5 and stiffness ratio kt/kn = 0.5. 

oriented contacts are proportional to those of eqn.(I7j), induced by the average 
strain and aligned with them. Hence we can write |32j : 



K{e) 
Me) 



cynKiO) 
at Am. 



The contributions from the particle displacements and rotations to A^ are 
characterized by analogous coefficients a" and af, such that 

ttt = a" + a^. 

Considerations analogous to those arisen from expression (??) allow one to 
write the effective shear modulus as 
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Figure 7: Magnitude of fluctuations in the relative displacements in terms 
of a^, at, a^ and a". Results from DEM simulations for three coordination 
numbers F: F = 3.5 (D), F = 4(A), F = 5 (Q) and various stiffness ratios 

5.3 Magnitude of observed and predicted fluctuations 

In the case of isotropic compression, Figure H] shows that the prediction of 
Pn is satisfactory. 

As regards shear, the dependence of a„, at, a" and a^ on coordination 
number and stiffness ratio kt/kn is shown in Figure [71 Negative numerical 
factors mean a relaxation with respect to the average strain assumption. 
The magnitude of the fluctuations decreases with increasing coordination 
number and, as expected, is more important when the system is undergoing 
shear deformation. Both the normal and the tangential component of the 
fluctuations relax the system with respect to the average strain assumption. 
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but along the tangential direction the relaxation is almost exclusively due 
to the particle rotation. On the contrary, the center displacements generally 
induce average relative displacements of the same sign as those due to the 
average strain. The predictions that include the fluctuations result in group- 



b -0 





Figure 8: Comparison between «„, at, af and af from DEM simulation (solid 
line), IPF (dashed) and PF-theory (dashed-dotted). Results for F = 3.5 and 
various stiffness rations kt/kn- 

averages analogous to those of Figure [71 The corresponding factors a„, at, 
a" and af are plotted in Figure [HI for coordination number F = 3.5. Their 
trend stays unchanged for the other coordination numbers considered. The 
estimates qualitatively reproduce the observed dependence of the fluctuations 
on kt/kn and on coordination number F. However, the relaxation they induce 
is underestimated, and the stiffening observed in the case of the particle 
displacement generally overestimated. The difference between estimate and 
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measured values at low coordination number in the case of shear emphasizes 
the necessity of considering larger correlation lengths if reliable predictions 
of the mechanical behavior of such systems are to be obtained. Given the 
remarkable improvement already obtained with respect to the average strain 
assumption, it might be sufficient to incorporate the fluctuations of the first 
neighbors. 

6 Discussion and perspectives 

This study has emphasized the insufficiency of the average strain assumption 
in predicting the elastic moduli of granular media, especially in presence of 
shear, as particles undergo important additional displacements and rotations 
in order to attain equilibrium. This particularly holds in the case of coordi- 
nation numbers smaller than that corresponding to the onset of iso-staticity 
for frictionless assemblies, whose occurrence seems relevant to practical pur- 
poses. 

We have shown that the role of the fluctuations is clearly phrased in 
terms of the relative displacements they induce at the contact points, coher- 
ently with the constitutive law for the contact forces. The relative displace- 
ments due to the fluctuations are highly correlated with contact orientation. 
Along the normal and the tangential direction, respectively, their average 
over equally oriented contacts is proportional to the relative displacements 
aligned with them and due to the average strain. Such a proportionality is 
expressed by the numerical factors /5„ for the case of isotropic compression, 
and an, «" and af in the case of shear. Their numerical value allows one to 
interpret the role of the different kinematic ingredients to the relaxation ob- 
served at the macroscopic scale. We have found that the normal component 
of the center displacements and the rotations always oppose the effect of the 
average strain. On the contrary, the tangential component of the center dis- 
placements generally stiffens the assembly, with the exception of large ratios 
kt/kn at large coordination number. 

We have analyzed two approaches that include the fluctuations from the 
average strain into the prediction of the elastic moduli, namely the "particle- 
fluctuation" and the "pair-fluctuation" method. They both determine the 
fluctuations by considering the problem of equilibrium of a small subassem- 
bly. That is, they are based on the assumption that the fluctuations organize 
with short correlation length, in the order of three or four diameters. The 
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comparison with the DEM simulations shows a significant improvement in 
predicting the moduh with respect to the average strain assumption. The de- 
formation mechanisms are quahtatively captured, even at small coordination 
numbers, together with their dependence on the ratio between tangential and 
normal stiffness and on coordination number, thus proving the correctness of 
the approach. However, the models do not capture with sufficient accuracy 
the fluctuations at low coordination numbers when a shear loading is applied 
(especially for the tangential relative displacements). This may be caused by 
the occurrence of larger correlation lengths than those assumed here. 

Future work will focus on identifying the correct correlation length. An 
other open issue is the mechanical behavior at ratios of tangential to normal 
stiffness larger than one, that easily occur at the contact between cylinders 
far from the onset of sliding. The independence of the a's on the contact 
orientation assesses that the fluctuations that determine the macroscopic 
behavior originate in the average geometry of the assembly. Therefore, ana- 
lytical approaches can capture them. As the average geometry is disordered, 
its representation has to be resolved in statistical terms, as done in the an- 
alytical solution in [IB]. Even though the issue has not been dealt with in 
this manuscript, we anticipate that the cited analytical prediction gives less 
accurate results than the numerical implementation of the corresponding ap- 
proach. This emphasizes the sensitivity of the modeling to the details of the 
statistical representation. This representation is still a critical issue whose 
improvement we will pursue. 
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